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ABSTRACT 

We present new equilibrium component distribution functions that depend on three 
analytic integrals in a Stackel potential, and that can be used to model stellar discs 
of galaxies. These components are generalizations of two-integral ones and can thus 
provide thin discs in the two-integral approximation. Their most important properties 
■ are the partly analytical expression for their moments, the disc-like features of their 

configuration space densities (exponential decline in the galactic plane and finite extent 
in the vertical direction) and the anisotropy of their velocity dispersions. We further 
^ ' show that a linear combination of such components can fit a van der Kruit disc. 
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pH , It has been known for a long time that the classical two-integral equilibrium theory in axisymmetric geometry is not sufficient 
O ^ to adequately describe the stellar discs of galaxies. In accordance with Jeans theorem (Jeans 1915), the phase space distribution 
function of a stellar system in a steady state depends only on the isolating integrals of the motion; the binding energy E and 
the vertical component of the angular momentum L z are isolating integrals in a stationary and axisymmetric system. It is 
a fundamental property of all two- integral distribution functions F(E,L Z ) that the dispersion of the velocity in the radial 
| direction equals the dispersion in the vertical direction: we know that, for example, the disc of the Milky Way does not have 
that property (Binney & Merrifield 1998). Illustrations of other shortcomings of a two- integral model in a galactic context 
can be found in Durand, Dejonghe & Acker (1996). 

The introduction of a third integral of the motion helps to overcome these constraints: in that case, the velocity dispersions 
5_J ■ can be different in all the directions. Numerical experiments show that a third isolating integral seems to exist for most orbits 
d " in realistic galactic potentials (Ollongren 1962, Innanen & Papp 1977, Richstone 1982). This third integral can be taken into 
account numerically in the models by using extensions of Schwarzschild's (1979) orbit superposition technique (Cretton et 
al. 1999, Zhao 1999, Hafner et al. 2000). It is also possible to define an analytic third integral specific to particular orbital 
families (de Zeeuw, Evans & Schwarzschild 1996, Evans, Hafner & de Zeeuw 1997) or an approximate global third integral 
(Petrou 1983ab, Dehnen & Gerhard 1993), but we choose to construct models with an exact analytic third integral by using 
a Stackel potential (Stackel 1890, de Zeeuw 1985). 

It is not quite obvious to define suitable global distribution functions F(E, L z , I3) that depend on three exact analytic 
integrals and that can somewhat realistically represent our ideas of a real stellar disc. For example, Bienayme (1999) made 
three-integral extensions of the two-integral parametric distribution functions described in Bienayme & Sechaud (1997), but 
these ones were built to model the kinematics of neighbouring stars in the Milky Way only. Dejonghe & Laurent (1991) also 
defined the three-integral Abel distribution functions, but these ones could not provide very thin discs in the two-integral 
approximation. Robijn & de Zeeuw (1996) constructed three-integral distribution functions for oblate galaxy models, but they 
also had problems to recover the two-integral approximation. 

In this paper we continue the work of Batsleer & Dejonghe (1995), who constructed component distribution functions 
that are two-integral, but that can represent (very) thin discs when a judicious linear combination of them is chosen. We 
use these components as a basis for new component distribution functions that are three-integral, of which the Batsleer & 
Dejonghe components are a special case. 

In the next section, we outline some fundamentals of two-integral equilibrium systems and we show how to model discs 
with a finite extent in the vertical direction. In section 3, we present some general facts about Stackel potentials and we 
present new analytic three-integral distribution functions that can represent stellar discs. An analytical expression for the 
moments of these distribution functions is calculated in section 4. In the next section, we discuss their physical properties and 
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show their realistic disc-like character. Finally, in section 6, we show that these distribution functions can be used as basis 
functions in the modeling of a van der Kruit disc. For the conclusions, we refer to section 7. 



2 TWO-INTEGRAL FUNDAMENTALS 

We denote the gravitational potential in cylindrical coordinates (jjj, <f>, z) by V(vj,z) = —ij)(w,z) < 0. The two classical 
isolating integrals of the motion are the binding energy, E = ip(zu, z) — \v 2 , and the z-component of the angular momentum, 
L z = ZUVrf, — zu 2 <f>. 

When we use the term orbit, we do not consider the information contained in the phases of the orbital motion: an orbit 
is thus shorthand for an orbital density. Even with this definition, each pair (E, L z ) represents a family of orbits; to uniquely 
identify one particular orbit, we need the presence of a third effective integral of the motion (i.e. we need a triple (E, L z ,Is)). 

The axisymmetric nature of the system implies we can focus on the motion in a meridional plane (i.e. a plane with 
constant <j)). For a position (zuo,zq) in the meridional plane, the expressions for E and L z imply that all families of orbits 
that visit this position have isolating integrals of the motion (E,L Z ) that meet the requirement 

T 2 

E<iP(w ,Zo)- r-^, (1) 
since we know that 

vl + v\ > 0. (2) 

For given E and L z , this relation restricts possible motion for the corresponding family of orbits to a toroidal volume in 
configuration space. 

In (E, L^)-space, Eq. ([j]) defines the region in which the points correspond to families of orbits passing through (zuo,zo). 
The boundary line (equality in Eq. ([!])) contains orbits that reach the given position with zero velocity (g) in the meridional 
plane. Keeping z = zo fixed while allowing zu to vary then gives us a family of such boundary lines, of which we denote the 
envelope by 

E = S Z0 (L Z ), (3) 
with the parametric equations (Batsleer & Dejonghe 1995, Eq. 4 & 5) 

{L 2 
E = ^(zu,z )-^ (4) 
L\ = -^ 3 £(^, 2 «,). 

All points in integral space with E < S ZQ (L Z ) represent families of orbits that will pass through z — zo at a certain zu. Orbits 
for which E = S ZQ (L Z ) also do reach the height Zq, but can never go any higher. All points in integral space with E > S ZQ (L Z ) 
represent families of orbits that cannot reach z = Zq. S zo is thus the minimal binding energy of an orbit that cannot bring a 
star higher than zo above the galactic plane. 

For every height z\ > zo we find a similar curve E = S Z1 (L Z ), with S Z1 (L Z ) < S ZQ (L Z ) for every value of L z . Similarly, 
the envelope for the orbits that cannot go higher than z = is given by E = Sq(L z ), which gives us all circular orbits in the 
galactic plane. 

Orbits belonging to a disc with a maximum height zo are thus given by (E,L Z ) for which S ZQ (L Z ) < E < So(L z ) (the 
shaded area in Figure 1). Batsleer & Dejonghe (1995) constructed disc-like component distribution functions with a finite 
extent in vertical direction by setting them equal to zero for E < S ZQ (L Z ). In order to fully understand the components that 
we present here, that paper should probably be considered as preparatory reading. 



3 CONSTRUCTION OF THREE-INTEGRAL COMPONENTS 
3.1 Spheroidal coordinates and Stackel potentials 

We will work in spheroidal coordinates, since these coordinates allow a simple expression for our (axisymmetric) Stackel 
potential. Spheroidal coordinates are given by (A, <j>, v ), with A and v the roots for r of 

= 1 a < 7 < 0, (5) 



t + a t + 7 

and (zu, <f>, z) cylindrical coordinates. The parameters a and 7 are both constant and smaller than zero. 

A potential is of Stackel form, if there exists a spheroidal coordinate system (A, (f>, v), in which the potential can be written 

as 

v M =-m^m, (6) 

A — v 
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Figure 1. The shaded area is the area 
in the (E, L^)-plane for which the cor- 
responding families of orbits cannot go 
higher than zo above the galactic plane. 
The distribution function of a stellar disc 
with maximum height zo is null outside of 
this area. 



for an arbitrary function /(r) = (r + "/)G(r), t = X,v. The function — G(A) then represents the potential in the z = plane. 

For this kind of potential, the Hamilton-Jacobi equation is separable in spheroidal coordinates, and therefore the orbits 
admit three analytic isolating integrals of the motion. The third integral of galactic dynamics has the form 



More details can be found in de Zeeuw (1985) and in Dejonghe & de Zeeuw (1988). 



3.2 Modified Fricke components 

As mentioned before, we intend to create three-integral stellar distribution functions, for the construction of stellar discs: we 
want to achieve an exponential decline in the mass density for large radii, while we want to introduce a preference for (almost) 
circular orbits. 

It has been known for some time that two-integral models can describe very thin disc systems, with the restriction 
that both vertical and radial dispersions are equal (Jarvis & Freeman 1985). So we want to create three- integral distribution 
functions that can describe very thin discs in the two-integral approximation, unlike the Abel distribution functions (Dejonghe 
& Laurent 1991). 

The Fricke components (Fricke 1952) E V L^ favour that part of phase space where stars populate circular orbits, so they 
could be taken as a starting point. However, they cannot be used in their basic form to model discs with a finite extent because 
they populate orbits which can reach arbitrary large heights: therefore, we will take as a starting point the components defined 
in Batsleer & Dejonghe (1995, Eq. 19). 

In order to make the components depending on the third integral ^3, we introduce the factor (p + qE + rh\ + sl:j) s in 
which the parameter s (and 8) will be responsible for the three-integral character of the components. The coefficients p, q, r 
and s can, in the most general case, be functions of L z . 

This leads us towards a general three-integral disc component of the form 



h = 



~(L* +L 2 y ) + (7 -a) 



2 G(\)-G(y) 



(7) 



A- v 




(8) 



if 




(9) 
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The distribution function is identically zero in all other cases. The function f(L z ) is defined as 

= 1 + l- aL , (2L 2 z S Z0 (L z )f(S Z0 (L z )r e -^^. (10) 

The parameter a is the rotation parameter (the value of a influences only the odd moments of F, see section 4). If 
a — 0, there is no rotation for the component, and if a = +oo, F represents a maximum streaming component with no 
counter-rotating stars. 

The requested exponential decline in the mass density with large radii is controlled by the parameter 02 (and to some 
extent by the parameter ai, see section 5). 

The parameter r\ is responsible for the favouring of almost circular orbits, i.e. orbits with a binding energy E as close as 
possible to that of circular orbits in the galactic plane (see section 5). 

Furthermore, if we want to favour (almost) circular orbits, we have to suppose the distribution function to be an increasing 
function of E: this forces q to be positive. 

Since a large I 3 implies that the orbit can reach a large height above the galactic plane (see de Zeeuw 1985 for a complete 
analysis of the orbits in a Stackel potential), the orbits with small /3's have to be favoured in order to describe thin discs: 
this forces s to be negative. 

Other constraints (on p, q, r, s and 77) will be imposed in section 4, in order to enable the analytical calculations of the 
moments. 



4 MOMENTS 
4.1 Theory 

The moments of a distribution function F at the point (A, <f>, v) of a spheroidal coordinate system are defined as 

Ht,m,n(\v) = J J J F(E,L z ,I 3 )v l x v™v™dv x dv 4 ,dv„, (11) 

with v\, v<f,, v u the components of the velocity in the A, the <f> and the v direction of the spheroidal coordinate system, and I, 
m, n integers. 

The mass density, the mean velocity and the velocity dispersions of the stellar system represented by F can easily be expressed 
in terms of the moments ([n]) by 

p(A, v) = /xo,o,o(A, v) 
p{v</,)(\, v) = po,i,o(A, v) 
pa\{\v) = /ia,o,o(A,i/) 
P{ v 1){\ v ) = a»o,2,o(A,i/) 

pal(X, v) = ^0,0,2 (A, v) (12) 

To obtain the value of one of these moments, we have to integrate over the volume in velocity-space corresponding to all 
orbits that pass through the point (A, <j>, v). 

Since all three integrals of the motion are quadratic in v\ and u„, if I or n is an odd integer, the moment fj,i >m ,n is identically 
zero. If I and n are even integers, the moment fii im , n can be written as an integral computed in the integral space. We assume 
in sections 4.1 and 4.2 that a — (in that case, there is no rotation and if m is odd, the moment is zero) and that m is an 
even integer (the general case will easily be derived from this one in section 4.3). Under these assumptions, 

Rm,n = 2 ~ +1 i+n [%&v?dLl [ fdEdlJ ^~ S ^ Y (p + qE + rLl + aI 3 )'(lt- h^ih -Ia) 3 ^ (13) 
w{\ — v) t - J VM v J J \So-S Zo J 

where 7^" and are given by 

f I+(E,Ll) = (X + 7 )[G(X)-E]- ^L.Ll 

\ I-(E,Ll) = (v + 1 )[G(v)-E]-^ j Ll 

More details can be found in Dejonghe & de Zeeuw (1988). 

For our components given by Equation the integration surface in the (E, /3)-plane is defined by 

I3(E,L 2 Z ) <h <lt(E,Ll) 

S Z0 (L Z )<E<^~^ (15) 
p + qE + rL 2 z + sh > 0. 
The integration limits for the integral in L 2 Z will be determined in section 4.3. 



(14) 
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Figure 2. The integration area for the inner double integral in Equation (13) in the (E, /3)-plane (left panel, L z is fixed) and in the 
(x(E, y(E, 7"3))-plane (right panel, see Equation 19). 



4.2 (Partly) Analytical expression for the moments 

We want to reduce the triple integral (13) to a simple one by solving the innermost double integral analytically. The reader 
not interested in the mathematical details might step directly to section 5. 

For this double integral to be analytically solved, however, we will have to make use of those combinations of p, q, r and 
s (see section 4.3) for which the integration area is transformed into the triangle bounded by 

r h = i+{E,L 2 z ) 

I h=Is(E,Ll) (16) 



I p + qE + rhl + sh = 0, 

as shown in Figure []. In this situation, we can express the factor (E — S ZQ )/(So — S Zo ) as a linear combination of the other 
three factors in the integrandum (corresponding to the bounding lines of the integration surface): 

E ~ Sz ° = t (lt -h)+u{h-Is)+v{p + qE + rLl + sI 3 ). (17) 



So - S, 

We impose r\ to be an integer: then the double integral in the (J5, Ia)-plane transforms into a sum of integrals: 

53S(<)(''7 < ) ,,i tV ~ l ~ J uJ J J dEdh (p+i E + rL * + sI: ^ s+1 ( J 3 + - h) n+l -^- i - j {h - (is) 

For each integral in this summation, the integrandum consists of factors whose zero-points define the bounding lines for the 
integration surface in the (E, 73)-plane. These integrals can be solved analytically. 

We will be in this situation (for all the points (A, <j>, v) of configuration space, where the moments are calculated) whenever 
p + qE + rh\ + s/3 = does intersect the i?-axis for E > S Z0 {L Z ). 

In order to solve the integrals analytically, one uses the new integration variables x and y, defined by (for a fixed L 2 Z ) 

x(E,I 3 )=I+(E)-h nq) 
y(E,h) = h-Is(E). [W) 

The line in the (E, 73)-plane p + qE + rl? z +SI3 — becomes y = 2/ m ax(a;), and the root of j/ max (a;) = is x max (see Figure ^). 
To make a more compact notation possible, we first define the auxiliary functions g(r) and h(r), as 

g(r) = G(t) - ^^-r (20) 
2(r + a) 

and 

h(r) = a{T + i)-q. (21) 
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We then have 

= - jp + rLl - h{v)i>{\, v) + s(u + >y)G(u) - 2^(9 - ^ 2 )| (22) 

= ~~h(yj^ m (23) 
and 

ym^(x) = ^|(x max - x) (24) 
Solving the integral part of one of the terms in Eq. (IL8|) for y then yields 



I / max 



x^^—^dx / y^~ +i [h(\)(y - y max (x))] s+ > dy 



i) 



(A - v) s +< +1 J 



- (A _ (/p . i+T B(^+i + l,tf + i + l) y ^ + -- 3 |/ max ^)- +,+J+J+1 ^ (25) 
After solving for x (analogous to what we did for y), one obtains for the whole summation ( |li| ) 

r(u + i±i - i - j) r(=±i + i) r(<5 + % + i) 



r(5 + r> + i±2 +2 ) 

The coefficients t, u and v are calculated by equalizing term by term in equation (^) 
s^s^ = -i(A + 7) + u{y + 1 ) + vq 
= — t + u + vs 



(26) 



(27) 



§- =i (A + 7 ) G(A)-^fe) -«(" + 7) G(iz) 



2(i/+c«) 



+ vrL z + up 



We find for the coefficients 

t i ./ 

with 

/ ^(A, f) — Szo — , , 3?max , , , . . ; , . 

« = ,u = -ft(A)« - —,t =u +s(\-v)v (29) 

Now we have a one-dimensional numerical integration to perform, like for the Abel components (Dejonghe & Laurent 1991). 

4.3 The one-dimensional numerical integration 

The triple integral (13) is reduced to a simple one if we judiciously choose the parameters p, q, r and s. The double integral 
in the (E, 7a)-plane (for a fixed L z ) can be solved analytically whenever p + qE + rL 2 + s/3 = does intersect the i?-axis for 
E > S zo (Lz). It has to be the case for all the L z relevant in the integration. So if we take p = — S zo , q — 1 and r < 0, the double 
integral can be solved analytically because E — S zo + rL 2 + s/3 = intersects the _E-axis for the value E — S Z(t — rL 2 > S zo . 
In that case, the double condition (^) becomes a simple one because the condition E — S ZQ > is automatically verified when 
E - S ZQ + rL\ + sl 3 >0 since s < and I 3 > 0. 

We now have to determine the integration limits of the simple integral in L z . Since conditions (|l|) and Q imply that 

Szo (L z )<E<ip-^L, (30) 

the integration limits for the integral in L 2 Z are, in a first time, determined by the intersections of S zg (L z ) and the line 
E = ip — in the (E, L 2 )-p\ane (see Figure 3). 

Furthermore, in order to have a non-degenerate (i.e. not empty) triangle in Figure |^, we must have 

> 0. (31) 

Since the condition ip — 5-% — S Z() > is automatically verified when condition (^) is verified, the equality in ( ^l| ) fixes the 
minimal and maximal angular momentum to take into account in the integration. 
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Figure 3. The (E, L2)-plane. The integra- 
tions limits in L\, {L z ) 2 min and 
must be between the intersections of the 
curve E = S ZQ (L Z ) and the line E = 

ip — Condition ( |3l| ) then definitively 
fixes these limits. Eq and E\ are the min- 
imal and maximal energy taken into ac- 
count in the integration. 



Knowing that L 2 Z = ro 2 «5, the resulting expression for the moment (with a = and m even) becomes 



= 2~ 



T(S + V +i^+2) J {vl) 



(vir + —(s zo ) 



0+ai 



(V Y 

V^max / 



-J. X) X! ( i 



f] — 1 \ ri 

V 

J 



xt 



(-h(\))S-{-h(u))^-. ._„ ,_„ 

r(* + * + i) r(t, + i±i-i-j)r(2±i+j) (-fcM)^.^ 



r(* + i) 



r(r7 + ^)r(^±i) (-h(\)y 



(32) 



This integration has to be performed numerically. 

In the general case a ^ 0, the expression (B3) is still valid for the even moments. When m is odd, the integrandum has 
to be multiplied by a factor A: 



A = 



1 + e - £ra l'V 



(33) 



5 PHYSICAL PROPERTIES OF THE COMPONENTS 

In this section, we show the realistic disc-like character of our stellar distribution functions: their mass density has a finite 
extent in the vertical direction and an exponential decline in the galactic plane, they favour almost circular orbits and their 
velocity dispersions are different in the vertical and radial direction. By varying the parameters, we can give a wide range of 
shapes to the components. 

In order to illustrate the role of the different parameters, we calculate the moments of many component distribution 
functions with different values for the parameters. 

As galactic potential, we use one of the Stackel potentials described by Famaey & Dejonghe (2001), that are extensions 
of the ones described by Batsleer & Dejonghe (1994). In the implementation of the theory we choose p = —S zo , 5 = 1, and 
r = 0. 



5.1 The parameter zo 

This parameter was introduced in order to impose a maximum height above the galactic plane for the disc-like component 
(Figure 0): indeed, when E > S ZQ , an orbit cannot go higher than z — zo, and the distribution function (0) is null for 
E < S Z g\L z ). In order to model samples of stars belonging to populations with different characteristic heights above the 
galactic plane, we can use a set of components with different values for this parameter. 
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Figure 4. Contour plots of the mass density (i.e. the moment /io,0,o) i n a meridional plane, for components with the parameters 
(ai , c*2, P, 5, r;, zq, s, a) = (3, 3, 1, 1, 2, Zo, —0.5, —5), with zq equal to 4 kpc (left panel) and 2 kpc (right panel) respectively (note the very 
different scale for vj and z). The discs become thinner with smaller zq, while the mass density is zero above z = zq (and even so below 
z = 20 because s ^ 0). In this figure and in the following similar ones, every contour corresponds to a density that is a factor of 10 
smaller than the next lower contour. 
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Figure 5. This figure displays the decline of the logarithm of the galactic plane mass density (in arbitrary units and scaling) of different 
components for varying a±. A rising cri helps to produce mass close to the center. The other parameters have the same values as in 
Figure 4 (with zq equal to 2 kpc) except that 02 = 0. 
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Figure 6. Logarithm of the configuration space density (in the galactic plane and in arbitrary units) of different components for varying 
Q2- The other parameters have the same values as in Figure 4 (with zo equal to 2 kpc) except that a\ is adjusted to built-in a given 
exponential decline. We see that the components have an exponential decline in the galactic plane and that a.2 is roughly the reciprocal 
of the component's scale length. 



5.2 The parameter oi\ 

The parameter a 1 enters Eq. (10) as the exponent of S ZQ : so, for non-negative values of cti, the factor where it appears will 
behave as a declining function of L z , in the same way as S Z0 {L Z ) does, showing a steeper decline for larger qi (see Figure 
|H|). A large a\ thus results in a distribution function that favours a large fraction of bound orbits. When it is increasing, this 
parameter helps to produce mass close to the center. 

When a given exponential decline is requested, ct\ will be a function of the other parameters rather than a fixed parameter 
(see section 5.3). 



5.3 The parameter CV2 

The parameter Q2 occurs as exponent in the distribution function's exponential factor . Increasing values of this parameter 
will contribute to the mass distribution near the center (Figure^), like in the ot\ case (but exponentially). 

On the other hand, our potential is approximately Keplerian at very large radii: this implies that, in the galactic plane, 

Ll^vj (34) 

and that (Batsleer & Dejonghe 1995) 

S Z0 (L Z )~-1~1 (35) 

So, at very large radii, Q2 is the reciprocal of the component's scale length, if the contribution of the other factors to the 
mass density does not vary much with respect to zu (for very large zu). In practice, it is often desirable to use components 
for which an exponential decline and a given scale length (as determined by Q2) is already built-in between two radii (say 
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Position CKpc ) 

Figure 7. Contour plots of the configuration space density (the mass density) in a meridional plane, for components with the parameters 
(ai,oc2,P,6,Tl,Zo,8,a) = (3, 3, f3, 1, 2, 2, -0.5, -5), with (3 = (top left), (3 = 0.5 (top right), (3 = 4 (bottom left) and = 6 (bottom 
right). We see that the maximum number of stars moves away from the center and that the mass is more concentrated in configuration 
space for an increasing f3. 



vo\ and vj-z). In such cases, the parameter a± is adjusted in such a way that it corrects for the non-constant behaviour of the 
other factors at large radii, making the global contribution of all factors (except the one in 02) constant at w\ and ZU2- 



5.4 The parameter ft 

For this parameter, there are two distinct cases: f3 = and f3 > 0. In the first case, the density is maximum in the center and 
falls off smoothly. In the latter case, the density is null in the center since L z = for zu — 0. In order to model real stellar 
systems, we need components with f3 — to have some mass in the center. However, in a real galaxy, the maximum number of 
stars occurs in the intermediate region where the bulge meets the disc: this justifies the utilization of components with f3 > 
when modelling real stellar systems. We see the maximum density moving away from the center when (3 is rising (Figure [jj). 
We also see on Figure (M) that an increasing j3 will concentrate the mass in a smaller region of configuration space. 



5.5 The parameter r\ 

For a given L z , the largest value of the factor (E — Sz ) v is obtained when the binding energy E = So corresponds to the 
circular orbits in the galactic plane (see Figure 1). So, the parameter r\ is responsible for the favouring of almost circular 



orbits: a larger rj implies a larger contribution of almost circular orbits (Figure 5.5) and thus a mass density located closer to 
the plane. 



5.6 The parameter s 

Condition Q E > S zo — s/3 implies that for ^3 7^ and a strictly negative s, the orbits cannot reach the height zo above 
the galactic plane. We see on Figure ([]) that the height zo is reached only in the case s — 0. Furthermore, since a large 73 
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Figure 8. For I 3 = and L z = 0.1, this 
figure displays the value of a component 
distribution function (in arbitrary units) 
between E = S ZQ and E = So for rj = 1 
and rj = 8. The other parameters have the 
same values as in Figure 4 (with zq equal 
to 2 kpc). We see that the proportion of 
circular orbits (close to E = So) is much 
higher for rj = 8. 



s. 



0.14 



0.16 



0.1S 



0.2 



s 



E 



corresponds to an orbit that can reach a large height, the factor [E 



S ZQ + sls) s favours orbits that stay low. So, by setting 



s more negative, we confine the orbits closer to the galactic plane. 

A very important property of our components is the possibility of introducing a certain amount of anisotropy in the 
stellar disc: if we denote by a z the dispersion of the velocity in the direction perpendicular to the galactic plane, and by <7 ro 
the dispersion of the radial velocity in the galactic plane, then any nonzero s will produce a ratio less than 1 (Figure 
p"c| ). The ratio is closer to unity in the center than in the outer regions: this indicates the physically realistic character of our 
components. For s — 0, we find a z = cr ra since we are dealing with a two-integral component again. 

5.7 The parameter 5 

A large 8 has partly the same effects as a large rj: it favours circular orbits. Furthermore, a large S augments the effects of 
the negative s and forces the stars to stay close to the plane by favouring low /3-values. As we can see on Figure (JTTJ) , a 
component with a larger S has more stars in the galactic plane but shows a steeper decline with respect to z. 



6 MODELLING 

The component distribution functions described in this paper are very useful as basis functions in the method described 
by Dejonghe (1989), in order to model any observable quantities (spatial mass density, velocity dispersions, average radial 
velocities on a sky grid,...). As an illustration, we present the application of the method to fit a given spatial density po(w, z) 
(see Batsleer & Dejonghe 1995 for a similar application in the two- integral approximation). We look for a linear combination 
of our components 



A 

that fits po(w, z), with A = (ai, 0.2, j3, 8, rj, Zo, s, a) and ca the coefficients that are to be determined. 

In practice, to find this linear combination we must introduce a grid {vJi,Zi) in configuration space and minimize the 
quadratic function in ca: 



i L \ a / 

This minimization, together with the constraint that the distribution function must be positive in phase space, is a problem 
of quadratic programming (hereafter QP) described by Dejonghe (1989). 




(36) 




(37) 
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Position CKpc ) 

Figure 9. Contour plots of the configuration space density (the mass density) in a meridional plane, for components with the parameters 
(a\ , ct2, f), 8, rj, zq, s, a) = (3, 3, 1, 1, 2, 2, s, — 5), with s = (top left), s = —0.5 (top right), s = —1 (bottom left) and s = —4 (bottom 
right). The height zg is reached only in the case s = 0. The more negative s, the more the mass is concentrated near the galactic plane. 



Here, we choose to adopt for po a spatial density which closely resembles that of a real disc, i.e. a van der Kruit law, for 
which the vertical disribution is a good compromise between an exponential and an isothermal sheet (van der Kruit 1988). 

Po(vj, z) oc exp ) secn (^~) ( 38 ) 

In order to have a zero derivative with respect to zn on the rotation axis, we adopt a mass density that follows closely the van 
der Kruit law, without a cusp in the center (see also Batsleer & Dejonghe 1995): 

po{m ' z) = TT^jt exp (~/|) sech (^) ' (39) 

with Hr and h z denoting the horizontal and vertical scale factor, respectively. 

Since the moments po,o,o are dependent on the potential of the galaxy (including the dark matter), we have to choose 
a potential for the galaxy that contains the stellar disc we want to model. We adopt a Stackel potential with three mass 
components that produces a flat rotation curve and that therefore is a candidate potential for a disc galaxy (Famaey & 
Dejonghe 2001, see also Batsleer & Dejonghe 1994). 

The actual modelling follows the same strategy as followed by Batsleer & Dejonghe (1995), which we briefly repeat here 
for easy reference. The first step in the actual modelling consists in the selection of a subset of components out of the (infinite) 
set of possible components. This subset is chosen so that certain features, that we suppose to be present in the stellar disc, 
such as circular orbits, are included. For example, we expect the mass density corresponding to a component to have an 
exponential behaviour close to the mass density we want to model. The QP program first minimizes the function (^) for 
one component Fa and chooses the component of the initial subset that produces the lowest minimum for that function ( |37| ) . 
Then the program iterates, selecting and adding at each iteration the component which, together with the components already 
chosen in a previous run, produces the best fit. Once the minimum of the x 2 -variable does not change significantly any more 
with the addition of extra components, the program is halted because too low a value for \ 2 could imply that the QP program 
starts producing a distribution function featuring unnecessary oscillations. 
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Figure 10. This figure displays the ratio -2*- of several components (in the Plane) for varying s. The other parameters have the same 
value as in Figure 0. The dependence of the components on the third integral induces anisotropy. 

As an example, we model a modified van der Kruit disc with ha = 3kpc and h z — 0.25kpc. Batsleer & Dejonghe (1995) 
already showed that a linear combination of two-integral components (with s = and 5 = 0) could fit such a disc, but with 
Cro — a z . In order to model real anisotropic velocity data in the future, the dependence on the third integral will be needed. 
We show that, by choosing components with /3 = 0, 1, 3, 5, 7; ai = 1; c*2 = 0.15, 0.3, 2; zo = 1, 2,4; r) = 1, 5, 10; s = 0, —0.5, — 1; 
5 — 0.01, 1, 4 and a = in the initial subset, a fit with components featuring s/0 and 5^0 can be obtained too (see Figure 
12). 

The fit is obtained for a linear combination of 25 components at 231 configuration space points (206 degrees of freedom). 
If we assume relative errors of 6%, we obtain for our minimum x 2 = 246, and the probability that a value of \ 2 larger than 
246 should occur by chance is Q(246/2, 206/2) ~ 0.1 (Abramowitz & Stegun 1972), which makes the goodness-of-fit believable 
(Press et al. 1986). 

By using Stackel dynamics to model a galactic disc, we construct a completely explicit and analytic distribution function, 
with an explicit dependence on the third integral. Figure (13) displays the distribution function obtained by QP in function 
of E, for L z = 0.1 and for two values of ^3 (^3 = and ^3 = 0.05). For 73 = 0, the distribution function is non-zero if 
Sz < E < So (with zo — 4kpc); for Lj, > 0, instead, the maximum value of E is the one corresponding to infinitcsimally thin 
short axis tubes and is smaller than 5*o. We see on Figure (13) that the distribution function is decreasing with increasing 
73 (particularly near E = S ZQ ), and that it has some clumps. These clumps at 73 — 0.05 are not discontinuities since the 
distribution function is a linear combination of continuous components. 

Many different three-integral distribution functions correspond to a given spatial density, and there is no guarantee that 
they will yield realistic velocity dispersions. It is a major result of this paper to show that it is possible to find a linear 
combination of our components yielding realistic velocity dispersions. Figure (14) displays the ratio in the galactic plane: 
at the radius corresponding to the solar position in the Milky Way (7.5-8.5 kpc), the classical value of — ~ 0.4 is obtained. 
The local maximum in the ^x. curve is due to the individual shapes of the velocity dispersions curves (Figure 14). 
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Figure 11. Decline of the logarithm of the configuration space density (in arbitrary units) as a function of the height above the Galactic 
plane at to = 1 kpc for varying 8. A rising <5 implies a steeper decline. The other parameters have the same values as in Figure 4 (with 
zq equal to 2 kpc). 



7 CONCLUSIONS 

In this paper, we have constructed new analytic three-integral stellar distribution functions yielding cr ra 7^ a z : they are 
generalizations of two-integral ones that can describe thin discs with the restriction that a m = a z (Batsleer & Dejonghe 
1995). 

We first reduced the triple integral defining their moments to a simple one, like in the Abel case (Dejonghe & Laurent 
1991), by making some assumptions on the parameters. Then we looked for the effects of the different parameters and showed 
the disc-like (physically realistic) features of our distribution functions: they have a finite extent in vertical direction and 
an exponential decline in the galactic plane, while favouring almost circular orbits. A very important feature induced by 
the dependence on the third integral is their ability to introduce a certain amount of anisotropy, by varying the parameters 
responsible for this dependence (s and 5). 

We finally showed that a van der Kruit disc can be modelled by a linear combination of such distribution functions with 
an explicit dependence on the third integral and a realistic anisotropy in velocity dispersions. This implies that they are very 
promising tools to model real data with cr w 7^ <j z (Hipparcos data for example) by using the quadratic programming algorithm 
described by Dejonghe (1989). This will provide information on the dynamical state of tracer stars in the Milky Way (or on 
external galaxies). 
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Figure 12. Data for a van der Kruit disc (hn = 3kpc and h z = 0.25kpc) were generated in the region kpc < to < 10 kpc and 
kpc < z < 1 kpc. The initial subset of components was made of the components with /3 = 0,1,3,5,7; ot\ = 1; 02 = 0.15,0.3,2; 
zq = 1,2, 4; r) = 1,5,10; s = 0,-0.5,-1; S = 0.01,1,4; a = 0. Components with non-zero s and non-zero 5 were selected by the QP 
program. The crosses indicate the data for z = Opc, 200pc, 400pc, 600pc, SOOpc, lkpc (from top to bottom), with error bars. The solid 
lines correspond to the mass density of the components linear combination at these heights. The only region of our disc where the fit 
deviates a little from the data is the 2 kpc central region: our components are primarily intended to describe the outer regions rather 
than the central region of galaxies since most of them have a supermassive black hole at the center, sometimes associated with deviations 
from axisymmetry (bar). 

REFERENCES 

Abramowitz M., Stegun LA., 1972, Handbook of mathematical functions, Dover, New York 

Batsleer P., Dejonghe H., 1994, A&A, 287, 43 

Batsleer P., Dejonghe H., 1995, A&A, 294, 693 

Bienayme O., 1999, A&A, 341, 86 

Bienayme O., Sechaud N., 1997, A&A, 323, 781 

Binney J. J., Merrifield M.R., 1998, Galactic Astronomy, Princeton Univ. Press, Princeton 

Cretton N., de Zeeuw P.T., van der Marel R.P., Rix H.-W., 1999, ApJS, 124, 383 

Dehnen W.D., Gerhard O.E., 1993, MNRAS, 261, 311 

Dejonghe H., 1989, ApJ, 343, 113 

Dejonghe H., de Zeeuw P.T., 1988, ApJ, 333, 90 

Dejonghe H., Laurent D., 1991, MNRAS, 252, 606 

de Zeeuw P.T., 1985, MNRAS, 216, 273 

de Zeeuw P.T., Evans N.W., Schwarzschild M., 1996, MNRAS, 280, 903 



© 2001 RAS, MNRAS 000, §-|l| 



16 B.Famaey, K.Van Caelenberg & H.Dejonghe 




Figure 13. For L z = 0.1 and a fixed 1$ (left paneh/3 = 0, riht panekijj = 0.05), this figure displays the values of the distribution 
function (corresponding to the fit obtained in Figure 12) as a function of E (for the bound orbits). For ^3 = 0.05, the maximum value 
of E is the one corresponding to thin tube orbits and is smaller than So- 
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Figure 14. The left panel displays the ratio -2^- in the galactic plane for the fit obtained in Figure 12. The right panel gives the shapes 
of the individual velocity dispersions curve <r ro (solid line) and a z (dotted line) in the galactic plane. 
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